Aula 11 - Análise de variância e correlação

Métodos Quantitativos Aplicados à Ciência Política

Frederico Bertholini

Formalizando Hipóteses

Um caso comum - Esta análise está correta?

Comparação de médias

A resposta é não. (por quê?)

Em todo caso, ela ilustra uma situação bem comum na prática, onde se deseja comparar médias. No caso, deseja-se comparar as taxas médias de mortalidade em cidades onde as armas são proibidas ou liberadas. Deseja-se testar se a média de homicídios em cidades onde armas são liberadas é menor que a média de homicídios em cidades onde armas são proibidas. (Como você colheria dados para esse estudo?)

As hipóteses, portanto, são:

\(H_0:\mu_L=\mu_P\)

\(H_A:\mu_L<\mu_P\)

Lembrando a base de trabalho

summary(dfe %>% dplyr::select(-id))

     media           faltas         turma               idade      
 Min.   :40.00   Min.   : 0.00   Length:60          Min.   :18.00  
 1st Qu.:70.00   1st Qu.: 2.00   Class :character   1st Qu.:19.75  
 Median :73.75   Median : 4.00   Mode  :character   Median :22.00  
 Mean   :74.38   Mean   : 4.25                      Mean   :25.23  
 3rd Qu.:80.00   3rd Qu.: 6.00                      3rd Qu.:29.00  
 Max.   :95.00   Max.   :10.00                      Max.   :49.00  
       interess        tempocup                  escola       estcivil 
 Secundário:24   não tem   : 4   Tudo privada       :20   Casado  :17  
 Principal :34   até 2h    : 3   Maior parte privada:15   Solteiro:42  
 NA's      : 2   de 2h a 4h:11   Maior parte pública:18   NA's    : 1  
                 de 4h a 6h:42   Tudo pública       : 7                
                 + de 6h   : 0                                         
                                                                       
     nota1           nota2      
 Min.   :39.00   Min.   :41.00  
 1st Qu.:67.50   1st Qu.:72.00  
 Median :72.75   Median :76.75  
 Mean   :71.99   Mean   :76.76  
 3rd Qu.:78.00   3rd Qu.:82.25  
 Max.   :95.00   Max.   :97.00

Diferença entre médias (amostras não pareadas)

\(H_0:\text{A média de notas de casados e solteiros é igual}\) ou \(H_0:\mu_c-\mu_s=0\) ou \(H_0:\mu_c = \mu_s\)

\(H_1:\text{A média de notas de casados e solteiros é diferente}\) ou \(H_1:\mu_c-\mu_s \neq 0\) ou \(H_1:\mu_c \neq \mu_s\)

Variável dependente: Notas

Variável independente: Situação conjugal

O que eu quero testar? Se a situação conjugal faz diferença na nota.

É efeito? Não! (Pearl, 2020) Inferência vs. Causalidade

Como testar na prática? Distribuições:

dfe %>% drop_na(estcivil) %>% 
  ggplot(aes(fill=estcivil,x=media,color=estcivil,group=estcivil)) +
  geom_density(color=NA,alpha=.65) +  
  geom_vline(data=. %>% group_by(estcivil) %>% summarise(media=mean(media,na.rm = T)),
             size=2,aes(xintercept=media,color=estcivil)) + 
  guides(color="none") + theme_ipsum_mod

Como testar na prática? Vamos construir intervalos

dfe %>% drop_na(estcivil) %>% 
  ggplot(aes(fill=estcivil,x=media,color=estcivil,y=estcivil)) +
  stat_summary(fun=mean, geom="point") + 
  stat_summary(fun.data=mean_ci, geom="errorbar", width=0.2) +
  geom_vline(data=. %>% group_by(estcivil) %>% summarise(media=mean(media,na.rm = T)),
             size=2,aes(xintercept=media,color=estcivil)) + 
  theme_ipsum_mod +theme(legend.position = "none")

Como testar na prática? Teste-t

t_test_results <- dfe %>%
  t_test(formula = media ~ estcivil,
         order = c("Casado", "Solteiro"))
t_test_results

# A tibble: 1 × 7
  statistic  t_df p_value alternative estimate lower_ci upper_ci
      <dbl> <dbl>   <dbl> <chr>          <dbl>    <dbl>    <dbl>
1     -1.03  40.4   0.311 two.sided      -2.53    -7.52     2.46

Outro exercício

\(H_0:\mu_\text{3joan} = \mu_\text{3joad}\)

\(H_1:\mu_\text{3joan} \neq \mu_\text{3joad}\)

Variável dependente: Notas

Variável independente: Turma (apenas 3joan e 3joad)

Intervalos

dfe %>% dplyr::filter(turma %in% c("3joan","3joad")) %>% 
  ggplot(aes(fill=turma,x=media,color=turma,y=turma)) +
  stat_summary(fun=mean, geom="point") + 
  stat_summary(fun.data=mean_ci, geom="errorbar",width=0.2) +
  geom_vline(data=. %>% group_by(turma) %>% summarise(media=mean(media,na.rm = T)),
             size=2,aes(xintercept=media,color=turma)) + 
  theme_ipsum_mod +theme(legend.position = "none")

Como testar na prática? Teste-t

t_test_results <- dfe %>% 
  dplyr::filter(turma %in% c("3joan","3joad")) %>% 
  t_test(formula = media ~ turma)
t_test_results

# A tibble: 1 × 7
  statistic  t_df p_value alternative estimate lower_ci upper_ci
      <dbl> <dbl>   <dbl> <chr>          <dbl>    <dbl>    <dbl>
1      2.37  36.0  0.0232 two.sided       8.31     1.20     15.4

Olhando no infer graficamente

# calculate the observed statistic
media_turmas <- dfe %>% 
  dplyr::filter(turma %in% c("3joan","3joad")) %>% 
  specify(media ~ turma) %>%
  calculate(stat = "t", order = c("3joan","3joad"))

# generate the null distribution with the theoretical t
distribuicao_teorica <- dfe %>% 
  dplyr::filter(turma %in% c("3joan","3joad")) %>% 
  specify(media ~ turma) %>%
  hypothesize(null = "independence") %>%
  calculate(stat = "t", order = c("3joan","3joad"))

Visualizando

# visualize the randomization-based null distribution and test statistic!
distribuicao_teorica %>%
  visualize(method = "theoretical") + 
  shade_p_value(media_turmas,direction = "two-sided") + 
  labs(title = "Distribuição teórica",x="Estatística t",y="Densidade")

Usando ggpubr

library(ggpubr)
dfe %>% dplyr::filter(turma %in% c("3joan","3joad")) %>% 
  ggerrorplot(x = "turma", y = "media",color = "turma",
              position = position_dodge(0.5)) +
  stat_compare_means(aes(label = paste0(..p.signif..," ou p = ", ..p.format..)),
                     method = "t.test") + theme(legend.position = "right")

Tamanho do efeito

O d de Cohen pode ser usado como uma estatística de tamanho de efeito para um teste t de duas amostras.

É calculado como a diferença entre as médias de cada grupo, dividido pelo desvio padrão agrupado dos dados.

Um d de Cohen de 0,5 sugere que as médias diferem pela metade do desvio padrão dos dados. Um d de Cohen de 1,0 sugere que as médias diferem por um desvio padrão dos dados.

dfe %>% dplyr::filter(turma %in% c("3joan","3joad")) %>%
  effectsize::cohens_d("media","turma",data = .)

Cohen's d |       95% CI
------------------------
0.77      | [0.10, 1.42]

- Estimated using pooled SD.

Testando diferentes hipóteses

Amostras pareadas

Uma mesma medição em dois momentos no tempo para os mesmos indivíduos

\(H_0:\mu_{t2}-\mu_{t1}=0\) ou \(H_0:\mu_{\Delta}=0\)

\(H_A:\mu_{t2}-\mu_{t1}\neq0\) ou \(H_A:\mu_{\Delta}\neq0\)

Ex.:: Nota 1 e Nota 2

Com infer

dif_mu <- dfe %>%
  mutate(d=nota2-nota1) %>%
  specify(response = d) %>%
  hypothesize(null = "point", mu = 0) %>%
  calculate(stat = "t")

hipotetica <- dfe %>%
  mutate(d=nota2-nota1) %>%
  specify(response = d) %>%
  hypothesize(null = "point", mu = 0) %>%
  calculate(stat = "t")

hipotetica %>% visualize(method="theoretical") + 
  shade_p_value(dif_mu,direction="two-sided") +
  labs(title = "Distribuição teórica",x="Estatística t",y="Densidade")
dfe %>% mutate(d=nota2-nota1) %>%  t_test(d ~ NULL)

Visualizando

# A tibble: 1 × 7
  statistic  t_df  p_value alternative estimate lower_ci upper_ci
      <dbl> <dbl>    <dbl> <chr>          <dbl>    <dbl>    <dbl>
1      10.8    59 1.30e-15 two.sided       4.77     3.88     5.65

Por que ANOVA?

\[ H_{0}: \mu_{1}=\mu_{2}=\cdots=\mu_{k}, \quad H_{A}: \mu_{i} \neq \mu_{j} \text{ para pelo menos um par } i \text{ e } j \]

O que é ANOVA?

Variabilidade dentro dos grupos = Soma dos Quadrados Dentro (SQD) \[ S Q D=\sum_{j=1}^{c} \sum_{i=1}^{n_{j}}\left(X_{i j}-\bar{X}_{j}\right)^{2} \] Variabilidade entre grupos = Soma de Quadrados Entre (SQE)

\[ S Q E=\sum_{j=1}^{c} n_{j}\left(\bar{X}_{j}-\overline{\bar{X}}\right)^{2} \] Variabilidade total = Soma Total de Quadrados (STQ)

\[ S T Q=\sum_{j=1}^{c} \sum_{i=1}^{n_{j}}\left(X_{i j}-\overline{\bar{X}}\right)^{2} \]

ANOVA

\(\text{STQ} = \text{SQE} + \text{SQD}\)

Fração da variabilidade explicada pelo grupo = \(\frac{\text{SQE}}{\text{STQ}}\)

É possível que, na população, as médias dos grupos sejam iguais e, por acaso, as médias das amostras sejam diferentes.

Quanto maior a variabilidade entre grupos (SQE) e menor a variabilidade dentro dos grupos (SQD), mais evidências teremos que as médias são diferentes na população.

Princípio: Teste F: \(\frac{\text{Variância entre grupos}}{\text{Variância dentro dos grupos}}\)

\(F = \frac{\text{MQE}}{\text{MQD}}\)

Na prática

\(H_0:\text{A média de notas das turmas é igual}\) ou \(H_0:\mu_\text{3joad}=\mu_\text{3joan}=\mu_\text{5joan}\)

\(H_A:\text{A média de notas de pelo menos uma das turmas é diferente}\) ou \(H_A:\mu_\text{3joad} \neq \mu_\text{3joan} \neq \mu_\text{5joan}\)

Variável dependente: Notas

Variável independente: Turma

Função aov

ANOVAtest <- dfe %>% aov(.,formula = media ~ turma)
summary(ANOVAtest)

            Df Sum Sq Mean Sq F value Pr(>F)  
turma        2    703   351.5   4.116 0.0214 *
Residuals   57   4867    85.4                 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Teste de Tukey

Tt <- TukeyHSD(ANOVAtest)
Tt$turma %>% as.data.frame() %>% 
  rownames_to_column() %>%  mutate(rowname = gsub("\n","",rowname)) %>%
  knitr::kable(col.names = c("","Dif.","Lim inf","Lim sup","p-valor"),
               digits=3,format = "latex")

ANOVA com infer

observed_f_statistic <- dfe %>% 
  specify(media ~ turma) %>%
  calculate(stat = "F")

dfe %>% 
  specify(media ~ turma) %>%
  hypothesize(null = "independence") %>%
  visualize(method = "theoretical") + 
  shade_p_value(observed_f_statistic,
                direction = "greater")

dfe %>% 
  specify(media ~ turma) %>%
  hypothesize(null = "independence") %>%
  generate(reps = 100, type = "permute") %>%
  calculate(stat = "F") %>%
  get_p_value(obs_stat = observed_f_statistic,
              direction = "greater")

# A tibble: 1 × 1
  p_value
    <dbl>
1       0

Anova com ggpubr 1

dfe %>% 
 ggboxplot(x = "turma", y = "media",color = "turma", palette = "npg")+
 stat_compare_means(comparisons = 
                      list(c("3joad","3joan"),c("3joan","5joan"),c("3joad", "5joan")), 
                    label.y = c(110, 105, 100))+stat_compare_means(label.y = 120)

Anova com ggpubr 2

# Multiple pairwise test against a reference group
dfe %>% 
 ggboxplot(x = "turma", y = "media",color = "turma", palette = "npg")+
 stat_compare_means(method = "anova", label.y = 120)+ 
 stat_compare_means(aes(label = ..p.signif..),method = "t.test", ref.group = "3joad")

Two-way ANOVA (dois fatores)

summary(ANOVAtest2 <- dfe %>% aov(.,formula = media ~ turma + interess))

            Df Sum Sq Mean Sq F value  Pr(>F)   
turma        2   1016   508.1   7.521 0.00131 **
interess     1      2     1.6   0.024 0.87840   
Residuals   54   3648    67.6                   
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
2 observations deleted due to missingness

Two-way ANOVA (dois fatores)

library(agricolae); HSD.test(ANOVAtest2, trt = c("turma","interess"),console = T)


Study: ANOVAtest2 ~ c("turma", "interess")

HSD Test for media 

Mean Square Error:  67.5601 

turma:interess,  means

                    media       std  r       se  Min  Max    Q25   Q50    Q75
3joad:Principal  82.77778  9.052317  9 2.739832 70.0 95.0 80.000 85.00 85.000
3joad:Secundário 78.88889  7.817360  9 2.739832 65.0 90.0 75.000 80.00 85.000
3joan:Principal  68.95833 12.082027 12 2.372764 40.0 80.0 66.875 72.50 77.500
3joan:Secundário 74.16667  4.082483  6 3.355595 70.0 80.0 70.625 73.75 76.875
5joan:Principal  73.07692  4.466758 13 2.279678 67.5 82.5 70.000 72.50 75.000
5joan:Secundário 73.33333  7.071068  9 2.739832 62.5 85.0 70.000 70.00 77.500

Alpha: 0.05 ; DF Error: 54 
Critical Value of Studentized Range: 4.178265 

Groups according to probability of means differences and alpha level( 0.05 )

Treatments with the same letter are not significantly different.

                    media groups
3joad:Principal  82.77778      a
3joad:Secundário 78.88889     ab
3joan:Secundário 74.16667     ab
5joan:Secundário 73.33333     ab
5joan:Principal  73.07692     ab
3joan:Principal  68.95833      b

Homogeneidade

plot(ANOVAtest2, 1)

library(car);leveneTest(media ~ turma * interess, data = dfe)

Levene's Test for Homogeneity of Variance (center = median)
      Df F value Pr(>F)
group  5  0.8897 0.4948
      52

Normalidade

plot(ANOVAtest2,2); shapiro.test(x = residuals(object = ANOVAtest2))


    Shapiro-Wilk normality test

data:  residuals(object = ANOVAtest2)
W = 0.9331, p-value = 0.003273

Qui-quadrado

Vamos olhar as relações entre:

Estado civil e interesse na disciplina.
Interesse na disciplina e turma.

Teste de independência entre variáveis categóricas.

\(H_0:\text{Variáveis são independentes}\)

\(H_A:\text{Variáveis não são independentes}\)

Tabelas de contingência e Qui-quadrado

Atenção à sua variável de interesse

dfe %>% drop_na(estcivil,interess) %>%
  tabyl(estcivil,interess)

 estcivil Secundário Principal
   Casado          5        12
 Solteiro         18        22

dfe %>% drop_na(estcivil,interess) %>%
  tabyl(estcivil,interess) %>% 
  janitor::adorn_percentages("col") %>%
  janitor::adorn_pct_formatting()

 estcivil Secundário Principal
   Casado      21.7%     35.3%
 Solteiro      78.3%     64.7%

No infer

qui_quadrado <- dfe %>% drop_na(estcivil,interess) %>%
   mutate_if(is.factor,as.character) %>% 
  specify(interess ~ estcivil,success = "Principal") %>%
  calculate(stat = "Chisq")

teorica_qui_quadrado <- dfe %>% drop_na(estcivil,interess) %>%
  mutate_if(is.factor,as.character) %>% 
  specify(interess ~ estcivil,success = "Principal") %>%
  hypothesize(null = "independence")

Visualizando

teorica_qui_quadrado  %>%
  visualize(method = "theoretical") + 
  shade_p_value(qui_quadrado,
                direction = "greater")

O teste

dfe %>% drop_na(estcivil,interess) %>%
  mutate(estcivil=as.character(estcivil),interess=as.character(interess)) %>% 
  infer::chisq_test(interess ~ estcivil)

# A tibble: 1 × 3
  statistic chisq_df p_value
      <dbl>    <int>   <dbl>
1     0.644        1   0.422

Tabelas de contingência e Qui-quadrado

dfe %>% drop_na(turma,interess) %>%
  tabyl(turma,interess)

 turma Secundário Principal
 3joad          9         9
 3joan          6        12
 5joan          9        13

No infer

qui_quadrado <- dfe %>% drop_na(turma,interess) %>%
  mutate_if(is.factor,as.character) %>% 
  specify(interess ~ turma) %>%
  calculate(stat = "Chisq")

teorica_qui_quadrado <- dfe %>% drop_na(turma,interess) %>%
  mutate_if(is.factor,as.character) %>% 
  specify(interess ~ turma) %>%
  hypothesize(null = "independence")

Visualizando

teorica_qui_quadrado %>%
  visualize(method = "theoretical") + 
  shade_p_value(qui_quadrado,
                direction = "greater")

O teste

dfe %>% drop_na(turma,interess) %>%
  mutate_if(is.factor,as.character) %>% 
  chisq_test(turma ~ interess)

# A tibble: 1 × 3
  statistic chisq_df p_value
      <dbl>    <int>   <dbl>
1      1.03        2   0.596

dfe %>% drop_na(turma,interess) %>%
  mutate_if(is.factor,as.character) %>% 
  chisq_test(interess ~ turma)

# A tibble: 1 × 3
  statistic chisq_df p_value
      <dbl>    <int>   <dbl>
1      1.03        2   0.596

Correlação

Três conceitos para a mesma idéia:

Covariância
Correlação
Coeficiente de determinação (\(R^{2}\))

São grandes em valor absoluto se houver forte relação linear

Guess the correlation

Covariância e correlação

\[ \operatorname{Cov}(X, Y)=\frac{\sum_{i=1}^{N}\left(X_{i}-\mu_{x}\right)\left(Y_{i}-\mu_{y}\right)}{N} \]

Uma área positiva indica associação positiva entre as variáveis.

Mas como saber se é uma associação forte ou fraca?
Qual a unidade de medida da covariância?

Para eliminar a unidade de medida das variáveis, podemos usar a padronização z. Desta forma, obtemos o coeficiente de correlação, que é a covariância com variáveis padronizadas. Este coeficiente varia de -1 a 1.

\[ \operatorname{Corr}(X, Y)=\frac{\sum_{i=1}^{N}\left(\frac{X_{i}-\mu_{x}}{\sigma_{x}}\right)\left(\frac{Y_{i}-\mu_{y}}{\sigma_{y}}\right)}{N} \]

Faltas e Nota 2

dfe %>% ggplot(aes(x=faltas,y=nota2)) +
  geom_point() + geom_smooth(method="lm")

Com ggpubr

dfe %>% ggscatter(x = "faltas", y = "nota2",add = "reg.line",conf.int = TRUE)+
  stat_cor(method = "pearson",label.x=7)

Com ggpubr 2

dfe %>% ggscatter(x = "faltas", y = "nota2",add = "reg.line",conf.int = TRUE,
          color = "turma", palette = "jco",shape = "turma")+
  stat_cor(aes(color = turma), label.x = 7,label.y=c(110,105,100))

Correlação

Um bom Guia

Vendo correlações graficamente

dfe %>% select_if(is.numeric) %>% cor() %>%
  corrplot::corrplot(.,method="number",type="upper",diag=FALSE )

Vendo correlações graficamente

Para mais

Modelo Linear Simples

Modelo linear simples estimado por mínimos quadrados ordinários (MQO)..